/*
 *-----------------------------------------------------------------
        Program to generate a set of random
	points in the plane. These points then
	undergo a relaxation procedure.
	
	Marcelo Walter - UBC 97707921

	July 1994
	Jan  1995
	17 Feb 1995 - Finished version alpha
	11 Mar 1995 - Include forms interface
	19 Sep 1995 - Changed the directory where I'm
		          saving and loading files to
		          /imager/project/marcelow/giraffe/forms/images
	20 Sep 1995 - Included the option for saving the pattern
                  as a ps file
	27 Sep 1995 - Forms functions handled through callbacks

	29 Sep 1995 - Start changes to make generic source code, ie,
	              without graphics dependencies

	29 Jan 1996 - Finish merging the interactive version
	              with the non-interactive one. A compiler
                  directive	sets which version to compile.
                  The option GRAPHICS invokes all forms-related
                  calls



    Thompson Peter Lied

    27 Jun 2001 - included version according thesis
                  Marcelo Walter
 *-----------------------------------------------------------------
 */

/*
 *----------------------------------------------------------------
 *	 INCLUDE
 *----------------------------------------------------------------
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/*#include <sys/time.h>*/

#include "common.h"

#ifdef GRAPHICS
	#include <forms.h>
	#include "sgi/inter.h"
#endif

#ifdef ENABLE_CUDA
	#include <cuda_runtime_api.h>
#endif

/*
 *-----------------------------------------------------------------
 * 	EXTERN GLOBAL VARIABLES
 *-----------------------------------------------------------------
 */
extern CELL 	*headC; 			/* declared in list.c */
extern CELL 	*zedC;  			/* declared in list.c */

extern BUCKET 	*head[BUCKET_SIZE_X][BUCKET_SIZE_Y]; /* declared in bucket.c */
extern BUCKET 	*zed[BUCKET_SIZE_X][BUCKET_SIZE_Y];  /* declared in bucket.c */

extern HEAP 	*myheap; 			/* declared in heap.c */

/*
 *-----------------------------------------------------------------
 *	LOCAL GLOBAL VARIABLES
 *-----------------------------------------------------------------
 */
/* second way to store the different adhesion forces between
   different cell's types */
float adhes[HOW_MANY_TYPES][HOW_MANY_TYPES] = {
  {0.5, 0.0, 0.0, 0.0},
  {0.0, 0.5, 0.0, 0.0},
  {0.0, 0.0, 0.5, 0.0},
  {0.0, 0.0, 0.0, 0.5},
};

/* store probabilities for switching types */
float sw[HOW_MANY_TYPES][HOW_MANY_TYPES] = {
  {1.0, 0.0, 0.0, 0.0},
  {0.0, 1.0, 0.0, 0.0},
  {0.0, 0.0, 1.0, 0.0},
  {0.0, 0.0, 0.0, 1.0},
};

/* ortcell holds the information about splitting rates for
   cells and probabilities for a given type to occur.
   These probabilities are only used when there is no
   external input experiment file */
CTYPE	ortcell[HOW_MANY_TYPES];


/* Modificaded 21/06/2001 */
/* Weighting factor for the displacements */
float wd;

/* anisotropic weight */
float wa;

/* weigth for repulsive radius */
float wrp;

/* seed for the random number generator */
int seed;

/* max and min values for the cell position */
float min_x=0.0, max_x, min_y=0.0, max_y;

/* global radius */
double radius, nradius, rsquare;

/* orientation in degrees for the anisotropy */
float orientation=0.0;

/* orientation vector */
double Ox, Oy, OxSq, OySq;

/* number of iterations between two time steps */
int iterations;

/* number of buckets in X and Y */
int bx,by;

/* size of a bucket in x and y */
float sizeX, sizeY;

/* my list of voronoi points, ie, the points
which form the voronoi polygons */
voronoiType *myVoronoiList=NULL;

/* my list of faces of voronoi polygons */
faceType *myFacesList=NULL;

/* keeps how many vertices the voronoi diagram has */
int numberOfVoronoiVert;

/* number of cells. npoints is just an update value for ncells
ie, ncells always keep the value given from the interaction
with the user while npoints is update as new cells are added */
int npoints,ncells;

/* number of cells by type */
int nCellsType[HOW_MANY_TYPES];

/* times for simulation and update */
int total_time, update_time;

/* Anisotropy effects */
double AniX, AniY, AniCommon;

/* holds the "go to" type. Used when changing manually
   the type of a given cell */
int gtType;

/* write output as an internal format file; to a ps format file
   or both */
OUTFILEMODE outmode;

/* flag variables. Specify whether or not the PS output
   file will have borders and/or be filled */
int borders, colorPS;

/* output ps file name */
char outpsfilename[256];

/* output internal format file name */
char outfilename[256];

/* input mode: external file, interactive (i.e., the user
   selects which cells will differentiate), non-interactive
   (the system decides which cells differentiate */
INMODE inmode=INTERACTIVE;

#ifdef GRAPHICS
  /* flags for initialization and voronoi drawing */
  byte init, voronoi_filled;
  /* flag for painting the cells by age or not */
  byte aging;
#endif

/*
 *----------------------------------------------------------------------
 * 	MAIN
 *----------------------------------------------------------------------
 */
int main(int argc, char *argv[])
{
  char inp_par_file[128];	/* input parameter file */
  char inp_exp_file[128];	/* input experiment file */
  char *s;			/* for parsing command-line parameters */
  byte par_flag=FALSE;		/* flag for a parameter file */
  int i;
  /*struct timeval pt;
  struct timezone zpt;*/
  /* monitor time for simulation */
  int initTimeSimulation, finalTimeSimulation, totalTimeSimulation;

#ifdef GRAPHICS
  FL_OBJECT *obj;
#endif
  
#ifdef ENABLE_CUDA
	cudaError_t cudaErr;
	int cudaDeviceCount;
	
	cudaErr = cudaGetDeviceCount( &cudaDeviceCount );
	
	if( cudaErr != cudaSuccess ) {
		printf("No CUDA capable device found.\n");
	} else
	{
		for (int device = 0; device < cudaDeviceCount; ++device) {
			cudaDeviceProp deviceProp;
			cudaGetDeviceProperties(&deviceProp, device);
			printf("Device %d : %s.\n", device, deviceProp.name);
		}
	}
	
#endif

  /* parse the command line options */
  for( i = 1; i < argc; i++ ){
    if( *argv[i] == '-'){
      s = argv[i] + 1;
      switch(*s)
	  {
		/* option to read a file with the cells positions */
		case 'e':
			i++;
			if( i >= argc ) break;
			strcpy(inp_exp_file,argv[i]);
			inmode = EXTERNAL;
			printf("Input exp file = %s\n",inp_exp_file);
			break;

		/* print the help or "usage" for the programa */
		case 'h':
			printf("Usage: cm -p <name of an input parameter file>\n");
			printf("Usage: cm -e <name of an input experiment file>\n");
			exit(1);
			break;

		/* option to read a file with the parameters */
		case 'p':
			i++;
			if( i >= argc ) break;
			strcpy(inp_par_file,argv[i]);
			par_flag = TRUE;
			printf("Input parameter file = %s\n",inp_par_file);
			break;

		default:
			printf( "Bad command-line parameter: %s\n", argv[i] );
			exit( 1 );
			break;
      }
    }
    else{
      printf( "Bad command-line parameter: %s\n", argv[i] );
      exit( 1 );
    }
  }
 
  /* initialize seed for random functions */
  /*gettimeofday(&pt,&zpt);
  initTimeSimulation = pt.tv_sec;
  seed = pt.tv_sec/pt.tv_usec;*/
  seed = -11;

  /* initializes data structures */
  init_data_struct();

  /* either read a parameter file or initialize
     the default values */
  if ( par_flag )
    inp_param(inp_par_file);
  else
    init_param();

  /* Init variables related with the repulsive radius */
  init_radius();
  
  /* init the orientation vector */
  compOrientation();
  
#ifdef GRAPHICS
  
  foreground();
  fl_init();
  create_the_forms();	/* created by fdesign */
  
  /* init the forms my routine */
  init_the_forms();
  
  fl_set_form_position( pafo, 10, -25);
  fl_show_form( pafo, FL_PLACE_POSITION, TRUE, "Control - CM Model (NEW)");
  
  /* init devices */
  init_devices();
  /* init graphics */
  init_graphics(max_x,max_y);
  
  /* I have to check here if there is an
     experiment file to be loaded. I have to
     call the forms "loading" procedure */
  if (inmode == EXTERNAL){
    handles_input_file(inp_exp_file);
    init_radius();
    print_info();
    assign_buckets();
    init=TRUE;
    voronoi_filled = TRUE;
    /* Compute the Voronoi */
    computeVoronoi();
    /* Draws a filled voronoi diagram */
    drawVoronoi(npoints,TRUE);
  }
  
  while(TRUE){
    obj = fl_check_forms();
    /* get mouse events out of the forms loop */
    if(obj == FL_EVENT) getinput();
  }
  
#else
  initing( inmode, inp_exp_file );

  aniso_effect();

  simulation();

  free_buckets_list( bx, by );

  /* Compute the Voronoi */
  fprintf( stderr, "Started computing voronoi...\n");
  computeVoronoi();
  fprintf( stderr, "Done!\n");

  /* writes the output file */
  fprintf( stderr, "Writing output files...\n");
  writeOutputFile( outmode, outpsfilename, outfilename );
  fprintf( stderr, "Done!\n");

  /* get final time for simulation */
  /*gettimeofday(&pt,&zpt);
  finalTimeSimulation = pt.tv_sec;*/
  totalTimeSimulation = finalTimeSimulation - initTimeSimulation;
  printf( "Total time %f minutes (ie, %d seconds)\n",
	  totalTimeSimulation / 60.0,
	  totalTimeSimulation );

  exiting();

#endif
	return 0;
}

/*
 *-------------------------------------------------------------------------
 *	This function initializes the data structures
 *------------------------------------------------------------------------
 */
void init_data_struct(void)
{
 
  /* init global list */
  create_global_list();

  /* init bucket lists */
  create_bucket_list();
  
  /* get space for the heap */
  myheap = HeapInit(HEAP_GROW_FACTOR);
}



/*
 *------------------------------------------------------------------------
 *  This function initializes parameters
 *  These are the default values, in case the user doesn't
 *  supply an input parameter file
 *  Adhesion types and probability of switching types
 *  are set to default values by initialization
 *  of the variable adhes, declared in main.c
 *------------------------------------------------------------------------
*/
void init_param(void)
{
	
  /* default max values for the cells positions */
  max_x = MAX_X;
  max_y = MAX_Y;

  iterations = NR;
  wd = WD;
  ncells = NCELLS;
  total_time = FINAL_TIME;
  update_time = TIME_UPDATE;
  
  wa = WA;
  wrp = WRP;
  outmode=INTERNAL; borders = TRUE; colorPS = FALSE;
  orientation = 0.0;

  /* inits original type of cell colors, division rates and initial
     probabilities */
  init_ortcell();

  /* set cell type C as the "go to" default type */
  gtType = C;


#ifdef GRAPHICS
  aging = FALSE;
  init=FALSE;
  voronoi_filled=TRUE;
#else
  inmode = NON_INTERACTIVE;
  strcpy(outpsfilename,"foo.ps");
  strcpy(outfilename,"foo.out");
#endif
  npoints=ncells;
  nCellsType[C] = nCellsType[D] = nCellsType[E] = nCellsType[F] = 0;
}

void init_radius(void)
{
  float area;

  /*printf("Started init_radius!\n");*/

  /* Compute area and radius of repulsion for relaxation */
  area = max_x * max_y / (float)ncells;
  radius = sqrt(area);

  radius *= wrp;
  nradius = radius;

  /* radius squared is used when computing the forces */
  rsquare = nradius * nradius;

  /* bx and by are the number of buckets in the x and y directions */
  bx = (int) max_x / ceil(nradius);
  by = (int) max_y / ceil(nradius);

  sizeX = max_x / (float) (bx);
  sizeY = max_y / (float) (by);

  /*printf("Finished init_radius!\n");*/
}

void compOrientation(void)
{
  double rad;

  rad = orientation * MINE_PI / 180.0;
  Ox = cos(rad);
  Oy = sin(rad);
  OxSq = Ox * Ox;
  OySq = Oy * Oy;
}


/*
 *-----------------------------------------------------------------
 *  Anisotropy effects. According Marcelow Thesis
 *----------------------------------------------------------------
 */
void aniso_effect (void)
{
  AniX = wa * Ox;
  AniY = wa * Oy;
  AniCommon = (1 - wa) * wd;
}


/*
 *------------------------------------------------------------------
 *	This function does the first setup for the simulation.
 *	Main tasks:
 *	- display initial configuration of cells
 *	- compute the first initial relaxation (longer than the others)
 *	- computes and draws the voronoi associated with the initial
 *	  set of cells
 *------------------------------------------------------------------
 */
void set_up(void)
{

  /* prints all points according to the buckets they belong to */
  /*print_buckets(bx,by);*/

  /* prints the global list */
  /*printf("List before relaxation!\n");
    print_list(headC,zedC);*/
  
  /* Calls the relaxation procedure. "iterations"
     is the desired number of relaxation steps  */
   /* relaxation(500); */
  
  /*printf("List after relaxation!\n");
    print_list(headC,zedC);*/
  
  /* Compute the Voronoi */
#ifdef GRAPHICS
  computeVoronoi();
  /* Draws a filled voronoi diagram */
  voronoi_filled = TRUE;
  drawVoronoi(npoints,voronoi_filled);
#endif

  /* Print Voronoi faces list */
  /*print_voronoi_faces_list(npoints);*/
}

/*
 *------------------------------------------------------------------------
 * I have to check that when I am including a new cell
 * the time for split can not be smaller than the 
 * present time
 *------------------------------------------------------------------------
 */
void handle_split(CELL *cell)
{
  BUCKET *bucket;
  CELL *new_cell;
  CELL_TYPE new_type;
  EVENT *motherCellEvent, *childCellEvent;
  int ix,iy;
  float rate;
  double nx,ny;
	
  /* increments the total number of cells */
  npoints++;
  
  /* This part here reduces the size of the radius as
     new cells are being added to the simulation. More
     cells mean that the average ideal area for a cell is smaller */

  /* compute the new radius */
  nradius = wrp * sqrt( max_x * max_y / (float) npoints);
  /* radius square is used when computing the forces */
  rsquare = nradius * nradius;

  /* I have to reorganize the bucket information
     as I start to include more cells */
  rate = (float) npoints / (float) ncells;
  if (rate >= 1.5){
    rearrangeBuckets();
    /*printf("rate = %f ncells = %d npoints = %d\n", rate, ncells, npoints);*/
    ncells = npoints;
  }
	
  /* increments the cell generation */
  cell->gen++;
	
  /* Computes the position for the child cell wx and wy are here
     for anisotropic effects (There was an old code here)*/

  /* Included 21/06/2001 */
  nx = cell->x + fran( -nradius*0.01, nradius*0.01 );
  ny = cell->y + fran( -nradius*0.01, nradius*0.01 );

  /* check if new position is out of boundaries
     and fix it */
  fixCoordinates( &nx, &ny );

  /* Check here if the new cell should be of a different
     type than its mother */
  new_type = sw_type( cell->ctype );

  /* creates child cell */
  new_cell = create_cell( nx, ny, cell->r, new_type, 0, FALSE);

  /* time for splitting of the child cell */
  new_cell->tsplit += cell->tsplit;

  /* update time for next split of mother cell */
  cell->tsplit += poidev( ortcell[cell->ctype].rateOfDivision, &seed );

  /* test for aging */
  new_cell->r = new_cell->tsplit;

  /* find the appropriate bucket */
  findBucket( new_cell->x, new_cell->y, &ix, &iy );
	
  /* insert cell information into the bucket */
  bucket = balloc();
  bucket->info = new_cell;
  insert_bucket( bucket, head[ix][iy], zed[ix][iy] );

  /* get space for next splitting events in the events queue */
  if (( motherCellEvent = (EVENT *) malloc(sizeof(EVENT))) == NULL )
    errorMsg("Problem with malloc in handle_split (main.c!");
  if (( childCellEvent = (EVENT *) malloc(sizeof(EVENT))) == NULL )
    errorMsg("Problem with malloc in handle_split (main.c)!");

  /* includes mother's next split event into the heap */
  motherCellEvent->whichEvent = SPLIT;
  motherCellEvent->eventTime  = cell->tsplit;
  motherCellEvent->info = cell;
  HeapInsert( myheap, motherCellEvent );
  /* includes child's next split event into the heap */
  childCellEvent->whichEvent = SPLIT;
  childCellEvent->eventTime =  new_cell->tsplit;
  childCellEvent->info = new_cell;
  HeapInsert( myheap, childCellEvent );

#ifdef GRAPHICS
  fl_set_counter_value(fd_number_cells,npoints);
  update_individual_number_of_cells();
#endif
}


/*
 *-----------------------------------------------------------
 * Switch a Cell, ie, makes the cell 'p' received as a type
 * 'new_type'
 *-----------------------------------------------------------
 */
void switch_cell(CELL *p, CELL_TYPE new_type)
{
  int i;
  CELL *s;

  /*printf("Found match with cell(%f,%f)\n",p->x,p->y);*/
  
  nCellsType[p->ctype]--;
  nCellsType[new_type]++;
  
  p->ctype = new_type;

  /*p->type.cell_type = new_type;
  
  p->type.color[R] = ortcell[new_type].color[R];
  p->type.color[G] = ortcell[new_type].color[G];
  p->type.color[B] = ortcell[new_type].color[B];
  p->type.rateOfDivision = ortcell[new_type].rateOfDivision;*/
  
  /*p->tsplit = poidev(p->type.rateOfDivision,&seed);*/
  p->tsplit = poidev( ortcell[new_type].rateOfDivision, &seed);

  /* test for aging */
  p->r = p->tsplit;
  
  p->gen = 0;
    
#ifdef GRAPHICS
  /* paint the cell on the screen */
  i=0;
  s = myFacesList[i].site;
  while( s != p) s = myFacesList[++i].site;
  draw_cell_filled(s,i);
  update_individual_number_of_cells();
#endif
}

/*
 *-----------------------------------------------------------------------
 * This function receives the information to create
 * a new cell. It allocates memory for the cell
 * and inserts the cell in the linked list.
 * INPUT: x and y -> cell's center position
 *	 r -> radius of repulsion for the cell
 *	 t -> cell type
 *	 g -> cell generation
 *	 b -> is a cell at the border (artificial) (TRUE or FALSE)
 *-----------------------------------------------------------------------
 */
CELL *create_cell(double x,double y,double r,CELL_TYPE t, byte g, byte b)
{
  CELL *p;

  /* b is a flag which indicates whether the cell is
     a "fake" cell at infinity just to draw a nicer
     voronoi. I'm not counting them */
  if(b == FALSE) nCellsType[t]++;
  p = palloc();
  p->x= x;
  p->y= y;
  p->fx = 0.0;
  p->fy = 0.0;
  p->gen = g;
  p->ctype = t;
  p->r= r;
  p->border = b;
  
  /*p->type.color[R] = ortcell[t].color[R];
  p->type.color[G] = ortcell[t].color[G];
  p->type.color[B] = ortcell[t].color[B];
  p->type.rateOfDivision = ortcell[t].rateOfDivision;*/

  p->tsplit = poidev( ortcell[t].rateOfDivision, &seed );
  
  /* test for aging. I'll be using the radius
     field to hold the time the cell was created */
  p->r = p->tsplit;
  
  /*printf("x = %f y = %f gen = %d\n",x,y,g);*/
  insert( p, headC, zedC );
  
  return( p );
}

void check_limits_for_coordinates(double *px, double *py)
{
  if (*px > max_x) *px = max_x;
  else if (*px < min_x) *px = min_x;
  
  if (*py > max_y) *py = max_y;
  else if (*py < min_y) *py = min_y;
}


/*
 *-----------------------------------------------------------------
 *	Writes the output as a postscript file
 *	"mode" == TRUE -> fill the polygons
 *	"mode" == FALSE -> do not fill the polygons
 *-----------------------------------------------------------------
 */
void print_the_polygons1(FILE *fp, byte mode)
{
   int number_of_vert,i,j,k;
   CELL *p;
   double d[XY];
   CELL_TYPE current_type=C;

   /* sets the color for the current type */
   /*  if (mode == TRUE) fprintf(fp,"0.8 setgray\n"); */ /* Modificado em 10/10/2001 */
   if (mode == TRUE) fprintf(fp,"w\n");
   else fprintf(fp,"0 setgray\n");
   
   for(i = 0; i < npoints; i++)
   {
      number_of_vert = myFacesList[i].faceSize;
      // myFacesList é um array alocado em getVoronoi (likmyco.c:95)
      p = myFacesList[i].site;
      
      if ( mode == TRUE ){
         if( current_type != p->ctype ){
            switch( p->ctype ){
               case C:
                  current_type = C;
                  if (colorPS) fprintf(fp,"w\n");
                  else fprintf(fp,"o\n");
                  /*fprintf(fp,"0.8 setgray\n");*/
                  break;
               case D:
                  current_type = D;
                  if (colorPS) fprintf(fp,"b\n");
                  else fprintf(fp,"m\n");
                  /*fprintf(fp,"0.6 setgray\n");*/
                  break;
               case E:
                  current_type = E;
                  if (colorPS) fprintf(fp,"u\n");
                  else fprintf(fp,"n\n");
                  /*fprintf(fp,"0.2 setgray\n");*/
                  break;
               case F:
                  current_type = F;
                  if (colorPS) fprintf(fp,"1 1 1 setrgbcolor\n");
                  else fprintf(fp,"0.4 setgray\n");
                  /*fprintf(fp,"0.4 setgray\n");*/
                  break;
            }
         }
      }
      
      /* first polygon point */
      k = myFacesList[i].vorPts[0];
      d[X] = myVoronoiList[k].x;
      d[Y] = myVoronoiList[k].y;
      
      check_limits_for_coordinates(&d[X], &d[Y]);
      fprintf(fp,"%.2f %.2f moveto\n",d[X],d[Y]);
      
      /* rest of polygons points */
      for(j = 1; j < number_of_vert; j++) {
         k = myFacesList[i].vorPts[j];
         d[X] = myVoronoiList[k].x;
         d[Y] = myVoronoiList[k].y;
         check_limits_for_coordinates(&d[X], &d[Y]);
         fprintf(fp,"%.2f %.2f lineto\n", d[X], d[Y]);
      }
      fprintf(fp,"closepath\n");
      if(mode == TRUE) fprintf(fp,"fill\n");
      else{
         fprintf(fp,"stroke\n");  /* command */
         /* draws cell's centers */
         fprintf(fp,"%f %f 1 d\n",p->x,p->y);
      }
   }
}	

void print_the_vertices(FILE *fp)
{
  int number_of_vert,i,j,k;
  double d[XY];
  int counter=0;

  /* print the x array */
  fprintf(fp,"/x [\n");

  for(i=0; i < numberOfVoronoiVert; i++){
    d[X] = myVoronoiList[i].x;
    d[Y] = myVoronoiList[i].y;
    check_limits_for_coordinates(&d[X], &d[Y]);
    fprintf(fp,"%.2f ", d[X]);
    counter++;
    if (counter == 10){
      counter=0;
      fprintf(fp,"\n");
    }
  }

  fprintf(fp,"] def\n\n");

  counter=0;
  /* print the y array */
  fprintf(fp,"/y [\n");

  for(i=0; i < numberOfVoronoiVert; i++){
    d[X] = myVoronoiList[i].x;
    d[Y] = myVoronoiList[i].y;
    check_limits_for_coordinates(&d[X], &d[Y]);
    fprintf(fp,"%.2f ", d[Y]);
    counter++;
    if (counter == 10){
      counter=0;
      fprintf(fp,"\n");
    }
  }
  fprintf(fp,"] def\n\n");
}

/*
 *-----------------------------------------------------------------
 *	Writes the output as a postscript file
 *	"mode" == TRUE -> fill the polygons
 *	"mode" == FALSE -> do not fill the polygons
 *-----------------------------------------------------------------
 */
void print_the_polygons(FILE *fp)
{
  int number_of_vert,i,j;
  CELL_TYPE current_type=C;
  CELL *p;

  /* sets the color for the current type */
  if (!borders) fprintf(fp,"/border {} def\n");

  for(i=0; i < npoints; i++){
    number_of_vert = myFacesList[i].faceSize;
    p = myFacesList[i].site;
    /* set the appropriate color */
      if( current_type != p->ctype ){
	switch( p->ctype ){
	case C:
	  current_type = C;
	  if (colorPS) fprintf(fp,"w\n");
	  else fprintf(fp,"o\n");
	  break;
	case D:
	  current_type = D;
	  if (colorPS) fprintf(fp,"b\n");
	  else fprintf(fp,"m\n");
	  break;
	case E:
	  current_type = E;
	  if (colorPS) fprintf(fp,"u\n");
	  else fprintf(fp,"n\n");
	  break;
	case F:
	  current_type = F;
	  if (colorPS) fprintf(fp,"1 1 1 setrgbcolor\n");
	  else fprintf(fp,"0.4 setgray\n");
	  break;
	}
      }

    for(j=0; j < number_of_vert; j++){
      fprintf(fp,"%d ",myFacesList[i].vorPts[j]);
    }
    switch(number_of_vert){
    case 6:
      fprintf(fp,"h\n");
      break;
    case 7:
      fprintf(fp,"i\n");
      break;
    default:
      fprintf(fp,"%d f\n", number_of_vert);
      break;
    }
    /* draws cell's centers */
    if(borders) fprintf(fp,"%.2f %.2f 1 d\n",p->x,p->y);
  }
}



/*
 *-----------------------------------------------------------------
 * handles_output_file
 *
 * There are two possible formats to save a file: internal
 * format and a postscript format. 
 * In the internal format we save all cell's positions
 * and types, together with the "state" of the system
 * when the simulation was saved.
 *
 * INPUT: pointer to a file name
 *	  mode (either "internal" or "ps")
 *-----------------------------------------------------------------
 */
void handles_output_file(char *out_file, OUTFILEMODE mode)
{
  FILE *fp;
  CELL *p;
  int i,j,k;
  int number_of_vert;
  double d[XY];
  char par_file[256];
		
  if( ( fp = fopen( out_file, "w" )) == NULL)
      errorMsg("Could not open output file to write!");
  
  if ( mode == INTERNAL ){ /* output file is in internal format */
    
    strcpy( par_file, out_file );
    strcat( par_file, ".par" );
    writeParamFile( par_file );
    
    fprintf( fp, "%d\n", npoints );
    p = headC->next;
    while( p != zedC){
      fprintf( fp,"%d %lf %lf\n", p->ctype, p->x, p->y );
      p = p->next;
    }
    
  } else if (mode == PS) {	/* output file is postscript */

    /* write the ps file header */
    save_ps_header(fp, npoints);

    /*
     * Marcelo: there are 2 possibilities for outputting
     * postscript. The one I have been using uses the stack
     * mechanism and produces smaller files, but not all
     * printers and drivers can deal with this mechanism.
     * The other way is the explicitly definition of all
     * polygons. Now I am NOT using the stck mechanism
     *
     * To use the stack mechanism make the following
     * changes:
     * a) comment the 2 lines:
     * 	print_the_polygons1(fp, TRUE);
     * 	if(borders) print_the_polygons1(fp, FALSE);
     * b) uncomment the lines:
     * 	print_the_vertices(fp);
     *  print_the_polygons(fp)
     */

    print_the_polygons1(fp, TRUE);
    if(borders) print_the_polygons1(fp, FALSE);

    /* this prints the vertices arrays, x and y */
    /*print_the_vertices(fp);*/

    /*
     * This prints the faces. If we want to see the
     * borders between cells we call the same
     * function twice with different parameters
     */
    /*print_the_polygons(fp);*/

    fprintf(fp,"showpage\n");
    
  }
#ifdef GRAPHICS
  else { /* mode is RGB SGI */
    output_view_window_as_rgb (out_file);
  }
#endif
		
  fclose(fp);
  fprintf(stderr, "File %s written...\n",out_file);

}

void handles_input_file(char *inp_file)
{
  FILE *fp;
  double x, y;
  int i, type;
  double area;

  nCellsType[C] = nCellsType[D] = nCellsType[E] = nCellsType[F] = 0;
  
  /* open input file to read */
  if( (fp=fopen(inp_file,"rb")) == NULL)
      errorMsg("Could not open to read input cells file inside handle_input_file!");
  
  fscanf(fp,"%d\n",&ncells);
  for( i = 0; i < ncells; i++ ){
    fscanf(fp,"%d %lf %lf\n",&type,&x,&y);
    create_cell(x, y, 0.0, (CELL_TYPE) type, 0, FALSE);
  }
  fclose(fp);
  npoints=ncells;
  
#ifdef GRAPHICS
  fl_set_counter_value(fd_number_cells,npoints);
  update_individual_number_of_cells();
#endif

}


void makes_all_cells(void)
{
  int i;

  nCellsType[C] = nCellsType[D] = nCellsType[E] = nCellsType[F] = 0;
	
#ifdef GRAPHICS
  if(inmode == NON_INTERACTIVE){
    for(i=0; i<ncells; i++)
	create_cell(fran(min_x,max_x),fran(min_y,max_y),0.0,pick_type(),0,FALSE);
  }
  else{
    for(i=0; i<ncells; i++)
	create_cell(fran(min_x,max_x),fran(min_y,max_y),0.0,D,0,FALSE);
  }
  fl_set_counter_value(fd_number_cells,ncells);
  update_individual_number_of_cells();

#else
  /*
   * Makes "ncells" and insert them into the list
   */
  for( i = 0; i < ncells; i++ )
      create_cell(fran(min_x,max_x),fran(min_y,max_y),0.0,pick_type(),0,FALSE);

#endif
}


CELL_TYPE pick_type(void)
{
  float x;

  x = ran1(&seed);

  if ( x <= ortcell[C].prob) return(C);
  /*else if ( x <= (ortcell[C].prob + ortcell[E].prob )) return(E);
    else if ( x <= (ortcell[C].prob + ortcell[E].prob + ortcell[F].prob )) return(F);*/
  else return(D);
	
}

/*
 *------------------------------------------------------------------------
 * Calculate a type according with the probability of a given 'type'
 *------------------------------------------------------------------------
 */
CELL_TYPE sw_type(int type)
{
  float x;
  
  x = ran1(&seed);
  
  if ( x <= sw[type][C]) return(C);
  else if ( x <= (sw[type][C] + sw[type][D])) return(D);
  else if ( x <= (sw[type][C] + sw[type][D] + sw[type][E])) return(E);
  else return(F);
}

/*
 *------------------------------------------------------------------------
 *  This function initializes the information about the cell's types
 *------------------------------------------------------------------------
 */
void init_ortcell(void)
{
  /* init info for cell type C */
  ortcell[C].color[R] = COR_TYPE_C_R;
  ortcell[C].color[G] = COR_TYPE_C_G;
  ortcell[C].color[B] = COR_TYPE_C_B;
  ortcell[C].prob = 0.3;
  ortcell[C].rateOfDivision = 4;
  
  /* init info for cell type D */
  ortcell[D].color[R] = COR_TYPE_D_R;
  ortcell[D].color[G] = COR_TYPE_D_G;
  ortcell[D].color[B] = COR_TYPE_D_B;
  ortcell[D].prob = 0.7;
  ortcell[D].rateOfDivision = 24;
  
  /* init info for cell type E */
  ortcell[E].color[R] = COR_TYPE_E_R;
  ortcell[E].color[G] = COR_TYPE_E_G;
  ortcell[E].color[B] = COR_TYPE_E_B;
  ortcell[E].prob = 0.0;
  ortcell[E].rateOfDivision = 0;
  
  /* init info for cell type F */
  ortcell[F].color[R] = COR_TYPE_F_R;
  ortcell[F].color[G] = COR_TYPE_F_G;
  ortcell[F].color[B] = COR_TYPE_F_B;
  ortcell[F].prob = 0.0;
  ortcell[F].rateOfDivision = 0;
}


void inp_param(char *inp)
{
  FILE *fp;
  
  if( (fp=fopen(inp,"rb")) == NULL)
      errorMsg("Could not open input parameter file to read!");
	
  /* inits original type of cell colors.
     The division rates and initial probabilities are also
     setup by this procedure but are later overwritten
     by the parameters inputed below */
  init_ortcell();//Não precisa ser chamada aqui

  fscanf(fp,"%f %f\n",&max_x,&max_y);
  fscanf(fp,"%d\n",&iterations);
  fscanf(fp,"%f\n",&wrp);
  fscanf(fp,"%d\n",&ncells);
  fscanf(fp,"%d %d\n",&total_time,&update_time);
  fscanf(fp,"%f %f %f %f\n",
	 &(ortcell[C].prob),
	 &(ortcell[D].prob),
	 &(ortcell[E].prob),
	 &(ortcell[F].prob));
  fscanf(fp,"%d %d %d %d\n",
	 &(ortcell[C].rateOfDivision),
	 &(ortcell[D].rateOfDivision),
	 &(ortcell[E].rateOfDivision),
	 &(ortcell[F].rateOfDivision));
  fscanf(fp,"%f %f \n",&wd,&wa);
  fscanf(fp,"%f %f %f %f\n",&adhes[C][C],&adhes[C][D],&adhes[C][E],&adhes[C][F]);
  fscanf(fp,"%f %f %f %f\n",&adhes[D][C],&adhes[D][D],&adhes[D][E],&adhes[D][F]);
  fscanf(fp,"%f %f %f %f\n",&adhes[E][C],&adhes[E][D],&adhes[E][E],&adhes[E][F]);
  fscanf(fp,"%f %f %f %f\n",&adhes[F][C],&adhes[F][D],&adhes[F][E],&adhes[F][F]);
  fscanf(fp,"%f %f %f %f\n",&sw[C][C],&sw[C][D],&sw[C][E],&sw[C][F]);
  fscanf(fp,"%f %f %f %f\n",&sw[D][C],&sw[D][D],&sw[D][E],&sw[D][F]);
  fscanf(fp,"%f %f %f %f\n",&sw[E][C],&sw[E][D],&sw[E][E],&sw[E][F]);
  fscanf(fp,"%f %f %f %f\n",&sw[F][C],&sw[F][D],&sw[F][E],&sw[F][F]);
  fscanf(fp,"%hd %hd %hd\n",
	 &(ortcell[C].color[R]),
	 &(ortcell[C].color[G]),
	 &(ortcell[C].color[B]));
  fscanf(fp,"%hd %hd %hd\n",
	 &(ortcell[D].color[R]),
	 &(ortcell[D].color[G]),
	 &(ortcell[D].color[B]));
  fscanf(fp,"%hd %hd %hd\n",
	 &(ortcell[E].color[R]),
	 &(ortcell[E].color[G]),
	 &(ortcell[E].color[B]));
  fscanf(fp,"%hd %hd %hd\n",
	 &(ortcell[F].color[R]),
	 &(ortcell[F].color[G]),
	 &(ortcell[F].color[B]));
  fscanf(fp,"%d %d %d\n", &outmode, &borders, &colorPS);
  fscanf(fp,"%f\n",&orientation);
  fscanf(fp,"%s\n",&outpsfilename);
  fscanf(fp,"%s\n",&outfilename);
  fclose(fp);
  
  npoints=ncells;


#ifdef GRAPHICS
  aging = FALSE;
  init=FALSE;
  voronoi_filled=TRUE;
#endif
}


void writeParamFile(char *out_file)
{
  FILE *fp;
  char cmfile[256];
  char psfile[256];
  int i, length;

  length = strlen(out_file);
  for(i=0; i < length-4; i++) cmfile[i] = out_file[i];
  cmfile[i] = '\0';

  length = strlen(cmfile);
  for(i=0; i < length-3; i++) psfile[i] = cmfile[i];
  psfile[i] = '\0';
  strcat(psfile, "ps");
  
  if( (fp = fopen(out_file,"w")) == NULL)
      errorMsg("Could not open output parameter file to write!");
  
  fprintf(fp,"%f %f\n",max_x,max_y);
  fprintf(fp,"%d\n",iterations);
  fprintf(fp,"%f\n",wrp);
  fprintf(fp,"%d\n",npoints);
  fprintf(fp,"%d %d\n",total_time,update_time);
  fprintf(fp,"%f %f %f %f\n",
	  ortcell[C].prob,
	  ortcell[D].prob,
	  ortcell[E].prob,
	  ortcell[F].prob);
  fprintf(fp,"%d %d %d %d\n",
	  ortcell[C].rateOfDivision,
	  ortcell[D].rateOfDivision,
	  ortcell[E].rateOfDivision,
	  ortcell[F].rateOfDivision);
  fprintf(fp,"%f %f \n",wd,wa);
  fprintf(fp,"%f %f %f %f\n",adhes[C][C],adhes[C][D],adhes[C][E],adhes[C][F]);
  fprintf(fp,"%f %f %f %f\n",adhes[D][C],adhes[D][D],adhes[D][E],adhes[D][F]);
  fprintf(fp,"%f %f %f %f\n",adhes[E][C],adhes[E][D],adhes[E][E],adhes[E][F]);
  fprintf(fp,"%f %f %f %f\n",adhes[F][C],adhes[F][D],adhes[F][E],adhes[F][F]);
  fprintf(fp,"%f %f %f %f\n",sw[C][C],sw[C][D],sw[C][E],sw[C][F]);
  fprintf(fp,"%f %f %f %f\n",sw[D][C],sw[D][D],sw[D][E],sw[D][F]);
  fprintf(fp,"%f %f %f %f\n",sw[E][C],sw[E][D],sw[E][E],sw[E][F]);
  fprintf(fp,"%f %f %f %f\n",sw[F][C],sw[F][D],sw[F][E],sw[F][F]);
  fprintf(fp,"%d %d %d\n",ortcell[C].color[R],ortcell[C].color[G],ortcell[C].color[B]);
  fprintf(fp,"%d %d %d\n",ortcell[D].color[R],ortcell[D].color[G],ortcell[D].color[B]);
  fprintf(fp,"%d %d %d\n",ortcell[E].color[R],ortcell[E].color[G],ortcell[E].color[B]);
  fprintf(fp,"%d %d %d\n",ortcell[F].color[R],ortcell[F].color[G],ortcell[F].color[B]);
  fprintf(fp,"%d %d %d\n",outmode,borders,colorPS);
  fprintf(fp,"%f\n",orientation);
  fprintf(fp,"%s\n",psfile);
  fprintf(fp,"%s\n",cmfile);
  
  fclose(fp);
  fprintf(stderr, "File %s written...\n",out_file);
}


/*
 *------------------------------------------------------------------------
 *	Prints general useful information about the simulation
 *------------------------------------------------------------------------
 */
void print_info(void)
{
  printf("\n---------------------------------------------------------\n");
  printf("\t\tPARAMETERS FOR THIS SIMULATION (NEW QUEUEING SYSTEM!)\n");
  printf("---------------------------------------------------------\n");
  printf("MAX_X = %4.2f MAX_Y = %4.2f\n",max_x,max_y);
  printf("Niterations = %d Ncells = %d\n",iterations,ncells);
  printf("Weight for the repulsive radius = %3.3f Radius = %3.3f\n",wrp,radius);
  printf("Number of buckets in x = %d and y = %d\n",bx,by);
  printf("Size of buckets in x = %2.2f and y = %2.2f\n\n",sizeX,sizeY);
  printf("Final Time = %d Update Time = %d\n",total_time,update_time);
  printf("Weight for the relaxation force = %2.2f wa = %2.2f \n",wd,wa);
  printf("Probabilities C = %5.3f D = %5.3f E = %5.3f F = %5.3f\n",ortcell[C].prob,
     ortcell[D].prob,ortcell[E].prob,ortcell[F].prob);
  printf("Division C = %d D = %d E = %d F = %d\n",ortcell[C].rateOfDivision,
     ortcell[D].rateOfDivision,
     ortcell[E].rateOfDivision,
     ortcell[F].rateOfDivision);
  printf("Adhesion CC = %2.2f CD = %2.2f CE = %2.2f CF = %2.2f\n",
     adhes[C][C],adhes[C][D],adhes[C][E],adhes[C][F]);
  printf("Adhesion DC = %2.2f DD = %2.2f DE = %2.2f DF = %2.2f\n",
     adhes[D][C],adhes[D][D],adhes[D][E],adhes[D][F]);
  printf("Adhesion EC = %2.2f ED = %2.2f EE = %2.2f EF = %2.2f\n",
     adhes[E][C],adhes[E][D],adhes[E][E],adhes[E][F]);
  printf("Adhesion FC = %2.2f FD = %2.2f FE = %2.2f FF = %2.2f\n",
     adhes[F][C],adhes[F][D],adhes[F][E],adhes[F][F]);
  printf("switch[C][C] = %2.2f [C][D] = %2.2f [C][E] = %2.2f [C][F] = %2.2f\n",
     sw[C][C],sw[C][D],sw[C][E],sw[C][F]);
  printf("switch[D][C] = %2.2f [D][D] = %2.2f [D][E] = %2.2f [D][F] = %2.2f\n",
     sw[D][C],sw[D][D],sw[D][E],sw[D][F]);
  printf("switch[E][C] = %2.2f [E][D] = %2.2f [E][E] = %2.2f [E][F] = %2.2f\n",
     sw[E][C],sw[E][D],sw[E][E],sw[E][F]);
  printf("switch[F][C] = %2.2f [F][D] = %2.2f [F][E] = %2.2f [F][F] = %2.2f\n",
     sw[F][C],sw[F][D],sw[F][E],sw[F][F]);
  printf("Colors:\n");
  printf("C [%d] [%d] [%d]\n",ortcell[C].color[R],ortcell[C].color[G],ortcell[C].color[B]);
  printf("D [%d] [%d] [%d]\n",ortcell[D].color[R],ortcell[D].color[G],ortcell[D].color[B]);
  printf("E [%d] [%d] [%d]\n",ortcell[E].color[R],ortcell[E].color[G],ortcell[E].color[B]);
  printf("F [%d] [%d] [%d]\n",ortcell[F].color[R],ortcell[F].color[G],ortcell[F].color[B]);
  printf("Outmode = %d borders = %d colorPS = %d\n",
     outmode, borders, colorPS);
  printf("Ox = %f Oy = %f orientation = %f\n", Ox, Oy, orientation);

#ifndef GRAPHICS
  printf("Output files: PS = %s INTERNAL = %s\n", outpsfilename, outfilename);
#endif

}

/*
 *-------------------------------------------------------------------
 *	Build the firts round of relaxation events
 *-------------------------------------------------------------------
 */
void update_relax_events_queue(int init_time, int final_time)
{
  EVENT *a;
  int i,time;

  if ((a = (EVENT *) malloc(sizeof(EVENT)))== NULL)
      errorMsg("Problem with malloc in update_relax_events_queue!");
	
  for(time=init_time;time<final_time;time++){
    /* Put on the heap the relaxation events */
    for(i=0; i < iterations; i++){
      a->whichEvent = RELAX;
      a->eventTime = time + 1.0*i/(float)iterations;
      a->info = NULL;
      HeapInsert(myheap,a);
    }
  }
  free(a);
}

/*
 *-------------------------------------------------------------------
 *	Build the first heap, ie, the first round of cells division
 *	and update the cell's data structure with the next
 * 	time to split	
 *-------------------------------------------------------------------
 */
void update_split_events_queue(void)
{
  EVENT *event;
  int i;
  CELL *cell;

  if ((event = (EVENT *) malloc(sizeof(EVENT))) == NULL)
      errorMsg("Problem with malloc in update_split_events_queue!");
		
  cell = headC->next;
  while( cell != zedC )
  {
    /* includes split event in the heap */
    event->whichEvent = SPLIT;
    event->eventTime = cell->tsplit;
    event->info = cell;
    HeapInsert( myheap, event );
    /* updates next split time */
    cell->tsplit += poidev( ortcell[cell->ctype].rateOfDivision, &seed );
    cell = cell->next;
  }
  free(event);
}


/*
 *-----------------------------------------------------------
 *  Fix the coordinates of a cell. If the coordinates
 *  are outside the (max_x, max_y) limit,  brings
 *  them into this limit in a toroidal fashion
 *-----------------------------------------------------------
 */
void fixCoordinates(double *x,  double *y)
{
  if(*x < 0) *x += max_x;
  *x = fmod(*x, (double) max_x);
  if(*y < 0) *y += max_y;
  *y = fmod(*y, (double) max_y);
}
